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FOREWORD 

A  computational  study  has  been  conducted  of  the  plume  created  by  an  underwater  explosion. 
Calculations  have  been  performed  with  a  compressible  and  an  incompressible  method.  The  latter 
method  was  found  to  be  more  economical  in  treating  the  relatively  long  tenn  phenomena  assrxiated 
with  the  explosion  plume,  which  takes  several  seconds  to  fonn  and  decay.  Both  techniques  suggest 
that  the  formation  of  an  explosion  plume  can  be  divided  into  five  different  phases;  cavity  fonnation. 
cavity  collapse,  venting,  jetting,  and  rebound.  Unfortunately,  the  internal  plume  stnicturc  predicted 
by  calculation  cannot  be  verified  at  this  time  due  to  the  absence  of  experimental  data  desenbing  the 
density  distribution  within  the  plume. 

This  work  was  supported  jointly  by  the  Naval  Surface  \V:u-fare  Center  Independent  Research 
Fund  and  the  Office  of  Naval  Re.search  Mathematical  Science.s  Division  through  Dr.  Richard  Lau 
under  contract  N0001492\VX24138. 
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ABSTRACT 

A  computati''inal  study  has  been  conducted  of  the  plume  created  by  an  underwater  explosion 
Calculations  have  been  perfonned  with  a  compressible  and  an  incompressible  method.  The  latter 
method  was  found  to  be  more  economical  in  treating  the  relatively  long  tenn  phenomena  associated 
with  the  explosion  plume,  which  takes  ses'eral  seconds  to  fonn  and  decay.  Both  techniques  suggest 
that  the  formation  of  an  explosion  plume  can  be  divided  into  five  different  phases:  cavity  fonnation. 
cavity  collapse,  venting,  jetting,  and  rebound.  Unfonunaiely,  the  internal  plume  structure  predicted 
by  calculation  cannot  be  verified  at  this  time  due  to  the  absence  of  experimental  data  descnbing  the 
density  distribution  within  the  plume. 
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INTKODL’CTION 


A  large  body  of  infonnation  exists  concerning  plumes  generated  by  underwater  explosions.  A 
review  of  much  of  this  information  can  be  found  in  References  !  through  3,  which  describes  the 
general  shape  of  the  resulting  plume  (its  height  and  width).  However,  available  data  does  not 
address  the  issues  o'"  plume  structure  "^d  the  density  distribution.  The  objective  of  this  report  is  to 
study  these  questions  using  numerical  simulation  techniques.  Aside  from  being  of  genera!  scientific 
interest,  the  results  from  this  study  can  be  used  to  evaluate  several  defense  concepts  which  utilize 
plumes  generated  by  explosives  for  ship  defense."' 

As  descriled  by  Cole’  and  Young,*  the  underwater  explosion  can  be  divided  into  two  distinct 
phases:  a  compressible  and  an  incompressible  one.  The  compressible  phase  occurs  early  after  the 
initiation  of  the  explosion  and  is  marked  by  the  outward  propagation  of  a  shock  uave  from  the 
c.xplo.sian.  This  shock  wave  reflects  off  the  water  surface,  leaving  an  upwclling  of  cavitaied  water 
behind  which  initiates  the  plume.  The  siiock,  which  travels  on  the  order  of  15(K)  nVsec.  leaves  the 
domain  of  interest  within  10  to  100  msecs.  For  the  remainder  of  the  problem,  where  the  time  scale 
is  on  the  order  of  seconds,  water  can  be  effectively  treateu  as  incompressible.  The  dominant  feature 
of  this  phase  is  the  pulsating  explosion  bubble  which  migrates  towards  the  free  surface.  This  bub¬ 
ble  starts  out  as  a  high  pressure  gas  pocket,  expanding  until  the  pressure  inside  of  it  is  less  than  that 
of  the  surrounding  water.  It  then  contracts  and.  under  the  influence  of  gravity  or  a  surface,  it  loses 
its  spherical  symmetry,  often  forming  a  torus  with  a  high  velocity  water  jet  flowing  through  its 
center.  Due  to  these  pulsations,  the  plume  generated  by  the  explosion  bubble  is  sensitive  to  the 
state  of  the  bubble  as  it  nears  the  surface. 

The  general  approach  has  been  to  model  the  underwater  explosion  [dienomena  using  two 
different  approaches;  a  compressible  and  an  incompressible  one.  The  compressible  approach  is 
based  on  the  algorithm  suggested  by  Collela,  et  al.^  Here,  the  Euler  equations  are  solved  for  both 
the  air  and  water  using  Simple  Line  Interface  Condition  (SLIC)  methodology^  to  capture  the  free 
surface.  This  technique  should  be  most  appropria'e  for  the  etu'ly  phases  of  the  problem.  An 
incompressible  method,  described  in  References  7  and  8.  has  been  employed  to  treat  later  times. 
This  approach  can  efficiently  and  accurately  predict  explosion  bubble  dynamics,  including  migration 
and  jetting.^ 

The  bulk  of  this  repon  concerns  the  incompressible  calculations,  which  can  be  efficiently 
applied  to  a  large  number  of  cases.  In  particular,  a  parametric  study  was  conducted  for  100  pounds 
of  TNT  exploded  at  depths  between  5  and  20  feet,  with  computations  carried  out  at  approximately 
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2-foot  intervals.  I'he  nunierical  aecurac)  o!  tlic  calcuhitiun'.  'Acre  assessed  b\  coniparitig  lesmts 
using  different  computational  grids.  I'hc  cunipressible  calculaiujiis  ucic  ap|?liiu  tu  fAs>  eases,  a 
300-pound  charge  at  depth  of  5  feet,  and  a  HKi-potmd  charge  in  shallow  uatcr  at  depth  of  6  5  icet. 
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CHAFFER  2 

DESCRIPTIONS  OF  THE  NUMERICAL  APPROACHES 

Two  different  algorithms  are  used  to  study  the  phenomena  described  in  the  previous  section. 
The  first  algorithm  is  based  on  a  model  in  which  the  water  is  assumed  to  be  incompressible,  while 
the  second  algorithm  allows  for  the  effects  of  compressibility. 


INCOMPRESSIBLE  MODEL 

The  “incompressible”  model  is  based  on  a  generalized  formulation  of  hydrodynamics  first  pro¬ 
posed  in  Reference  10.  This  formulation  uses  a  fi,\ed  spatial  domain  Q,  where  the  density  (3,  velo¬ 
city  u,  and  the  pressure  P  are  governed  by  the  mass  and  momentum  conservation  equations 

p,  V  {pu)  =  0,  (2-la) 

(pu),  +  V-(puu)  =  -pgk-VP,  (2-lb) 

subject  to  the  constraint 

P  ^  Po-  (2-lc) 

where  g  the  gravitational  constant,  -k  the  unit  vector  in  the  direction  of  the  gravitational  force,  and 

Po  is  the  constant  density  of  the  incompressible  liquid.  In  constructing  solutions  for  the  constrained 
system  (2-1),  it  is  assumed  that  liquid-on-bquid  collisions  behave  inelastically,  thereby  causing  a 
reduction  in  the  total  mechanical  energy  of  the  flow  field.  These  energy  losses  are  associated  with 
breakdowns  of  the  classical  theory  and  may  be  attributed  to  turbulence." 

The  density  field  divides  into  tw'o  time  varying  regions,  namely,  the  liquid  region  D(r)  = 
{x  e  Q:  p(x,f)=po},  and  the  nonliquid  region  where  p  <  py.  The  interfaces  separating  these  regions 
are  the  free  surfaces.  The  numerical  solution  "captures"  these  surfaces  as  slightly  smetu-ed  inter¬ 
faces.  The  nonliquid  regions  are  characterized  by  specifying  uniform  pressures  in  each  of  its  con¬ 
nected  subsets.  For  e,sample,  in  the  atmospheric  region  above  the  liquid  the  pressure  is  set  to  the 
ambient  air  pressure.  In  an  underwater  bubble  the  pressure  Pjj  may  be  determined  using  the  adia¬ 
batic  law  Pjj  =  CVf'^,  where  Vq  is  the  instantaneous  bubble  volume,  C  is  a  constant,  and  y  is  the 
ratio  of  specific  heats  of  the  bubble  gases. 

Assume  that  the  density  and  velocity,  p”,u”  at  time  step  n,  and  the  pressure  gradient  at  the 
previous  half  step,  are  known.  This  solution  is  evolved  using  the  following  three  step  tunc 

split  procedure. 
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Convection 

The  solution  is  first  advanced 

.u” ) — ^(p.u) 

by  “solving”  the  conserv'tuion  laws  (2-la,b)  without  including  the  VP  tcnn  on  the  nghi  liand  side 
of  (2- lb)  and  without  regard  to  the  constraint  (2-lc).  This  step  is  implemented  numcncaliy  using  a 
formally  second  order  Godunov-type  method,  which  uses  slope  limiting  in  space  and  explicit 
predictor-corrector  time  stepping. 

Redistribution  of  Density  and  Momenta 

Next,  the  density  and  momenta  arc  redistributed 

(p,u^-^(p,u) 

so  that  the  constraint  (2-lc)  is  satisfied  ,  the  global  conservation  of  mass  and  momenta  are  main¬ 
tained,  and  the  energy  is  nonincreasing.  The  density  is  redistributed  using  an  approximate  solution 
to  an  obstacle  problem.  This  solution  is  obtained  using  a  constrained  direction  preconditioned  con¬ 
jugate  gradient  method.  The  momenta  redistributions  are  determined  as  solutions  of  two  elliptic 
self-adjoint  problems.  Discretizations  of  these  problems  yield  diagonally  dominant  matrices  which 
are  efficiently  solved  by  a  diagonally  preconditioned  conjugate  gradient  method.  After  this  step 
p  =  p'’■^^  and  the  new  nonliquid  region  is  then  determined  along  with  the  pressure  in  each  of  its 
connected  subsets.  In  the  computational  space  the  new  liquid  region.  O'*"’  =  D(t"^^),  is  defined  to 
be  the  collection  of  grid  cells  C,  ^  such  that 

p;:;'>(l-ep)Po,  (2-2) 

where  is  the  density  in  cell  C,  A  typical  value  used  for  Cp  is  0.04.  At  this  stage  u  = 
in  the  new  non-liquid  region.  However,  u  is  not  consi.sient  with  (2- lb)  in  the  new  liquid  region. 

Pressure  Projection 

In  this  step  the  velocity  is  corrected 

u  -> 

using  the  gradient  of  the  new  pressure,  p The  pressure  P  =  is  the  solution  of 

tap  =  PqV  u  in 

where  x  is  the  time  step  This  equation  is  discretized  using  a  finite  element  method  with  bilinear 
elements,  and  the  resulting  linear  system  is  solved  using  an  incomplete  Cholesky  preconditioned 
conjugate  gradient  method.  The  new  velocity 


Po 

is  divergence  free  in  D”""*  and  thus  is  consistent  with  (2- lb). 
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A  more  complete  description  of  the  motivation  and  discretization  of  the  above  three  steps  are  dis¬ 
cussed  in  References  7  through  9. 


COMPRESSIBLE  MODEL 


The  compressible  model  uses  the  technique  outlined  in  Reference  5.  Here,  the  How  is 
described  by  the  2-D  Euler  Equations  in  cylindrical  coordinates; 


dQ  dU  dF  ^ 

dt  ^  dr  dz 


(2-3) 


where 


p 

pll 

,  U  = 

1 - 

s. 

pv 

piiv 

[pE. 

pH  J 

pu 

pu 

F  = 

puv 

o 

,  s  = 

1 

pir  +  p 

pv-  -t-p 

r 

puv 

pH 

piiH 

E  =  e  +  .5{ir  '■  v*)  ,  H  ~  h  +  .5{u^  +  v‘). 


(2-4) 


Here,  e,  h,  u,  v,  p,  and  p  are  the  energy,  enthalpy,  radial  velocity,  axial  velocity,  and  density  in 
single  property  cells.  In  mixed  cells,  these  are  mass  averaged  properties  for  the  cell.  Additional 
variables  are  also  computed  in  mixed  cells:  /  (volume  fraction),  p  and  E  for  each  substtincc. 

The  numerical  technique  used  to  solve  the  above  equations  is  based  on  a  time  splitting 
approach  which  divides  the  two-dimensional  problem  into  two  one-dimensional  problems  in  the  r 
and  z  directions: 

...  ...  dQ  dU  „  r 

(r  direction)  -z—  +  -z—  =  S  (2-5) 

dt  dr 

(z  direction)  -f  —  =  0.  (2-6) 

dt  dz 

The  complete  problem  was  simulated  by  solving  the  above  at  each  step  in  succession 

A  second  order  Godunov  scheme  is  used  to  solve  each  one-dimensional  problem.  This  method 
consists  of  a  predictor  and  corrector  step  and  is  described  in  detail  in  Reference  12,  In  regions  of 
the  flow  containing  a  single  species,  this  method  consists  of  the  following  steps; 

1.  Property  derivatives  are  computed  using  a  limited  second  order  differencing  procedure,  that 
reduces  derivative  values  in  non-smooth  regions  of  the  flow  field. 
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2.  The  predictor  step  is  based  on  a  method  of  characteristics  and  determines  properties  at  the  ceil 

edges  center  (i.e.,  t  +  This  produces  two  sets  of  property  values  at  each  cell  edge 

center,  one  from  each  of  the  adjacent  cells 

3.  The  corrector  step  applies  a  linear  Riemann  problem,  using  these  two  sets  of  values  as  inilia’ 
states,  to  determine  corrected  properties  at  the  cell  edge  centers. 

4.  Corrector  cell  wdge  center  properties  are  used  to  advance  Q  and  hence  the  cell  properties. 

In  multiple  fluid  regions,  the  algorithm  must  be  modified.  The  first  modification  consists  of 
applying  the  single  species  procedure  using  average  properties.  For  an  air-water  mixture,  the  aver¬ 
age  density  and  energy  follow  from  conservation: 

P  +  AvP... 

„  f aPa^a  t 

£  =  - ,  (2-7) 

P 

where 

=  volume  fraction  of  air , 
and 


/ 


=  volume  fraction  of  water 


Here,  the  subscripts  a  and  w  refer  to  air  and  water,  respectively  The  average  pressure  is  calculated 
by  defining  and  F,,,,  the  inverse  of  which  arc  a  measure  of  the  compressibility  of  air  and  water, 
respectively.  For  example,  for  substance  a. 


Noting  that  V''  is  proportional  to  v,  the  specific  volume,  assuming  an  iscntropic  variation  of  v  with 
p,  and  applying  the  definition  of  the  speed  of  sound, 

2  .‘pPa) 

ai  =  ^  =  -v^  —  .  (2-9.i 

,  Pa  J,v  ["'"aj  s 

yields 


^  a  Pu 


(2-10) 


The  definition  of  F  for  a  mixture  follows  bv  espansion: 

1  = f £_]  f ^  _f ic]  f 

F  bp  V  bp  V 

^  J  -  J  L  ^  J  <  ) 
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J  JL.  1 

i5;m 


6\',  r 


<'S\' 


V. 


i;. 

T 


^  61' , 


\p^^ 

1  6/j  !i  r 

K  ’  J  K 


6\' 


-I\.  - 


~f:^  I 


(2-!  n 


To  d^'tcnnir 
proui;v;c>  ihi 


to 
viciJi 


5p 


^  L 

r/"  r. 

c  tl.L'  u%crage  pre^su^e  for  a  mixture,  a  change  in  substance  volume  is  eoniputed  which 
t  same  pressure,  p .  in  each  component  Detimng  the  ciiange  m  pressure  of  a  component 
1/?^  - /3  ).  ctifurcin^  cs'rtscrs ation  of  vuUime  6'>\.  -r  61',^  ~  (h  and  using  hquation  (2-K) 


h .  t  p , 


P) 


0, 


t2-12f 


.  J\JK.  ; 

=n-p—  4--^-  . 
!_  la  1.  j 

Finally,  the  average  speed  of  sound  for  a  mixt.  re  is  deSined  by 

‘1'  /h-C- 


(2-13) 


I'he  second  mudiheation  of  the  akorithm 


fur  mixtures  consists  of  the  following  steps: 


1. 


A 


Update  the  volum.e  fr:sct,on  of  air  and  'vat-‘r  in  each  cell.  T'/iis  is  accomplished  in  three  stages; 

a.  Detc.nr.itiC  the  distributiost  of  tiic  air  and  water  volume  fractions  in  each  cell  using  SL!C 
algorithm.  'I'his  method  examines  the  contents  of  neighboring  cel!''  and  sets  the  air-water 
inierf.icc  to  one  of  the  options  sliov.  n  in  Figure  2-1 

b  Use  the  cell  edge  velocity  detennined  hy  the  linear,  average  property  Rictnann  solution  to 
compute  the  volume  of  each  substance  crossing  a  cel!  edge, 

c.  Update  the  volume  fractions  in  each  cell  using  the  fluxes  determined  at  stage  b.  The  sum 
of  the  new  fractions  will  not  necessarily  equal  the  cel'  volume.  Adjustments  are  made  by 
compressing  or  expanding  the  cells’  contents  assuming  a  reiativ  '.•ompressibility  for  air 
end  water  as  defined  by  l.'r^,  and  IfU^,. ,  The  resulting  correction  to  f  is 


,  r  1 

/« 

1  -1-  — 1 

V'  -  V'! 

1 

L  J 

(2-14) 


where  quantities  with  a  tilde  are  the  uncorrccted  values. 

Compute  the  air  and  water  mass,  momentiun  and  energy  flux  into  each  ceil  This  is  accom¬ 
plished  by  applying  first  order  upwind  differencing.  Using  the  volume  frtictions  eoniputed  in 
die  previous  step,  the  density  and  energy  of  each  substance  can  be  determined, 

Detenrune  an  average  cell  pressure  using  Equation  (2-12).  Adjust  the  volume  fraction  of  each 
substance  to  make  its  pressure  ctjual  to  the  average  pre.ssurc  using 
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1/^u  -P)  i 
Vp  J 


4.  Alter  Pa  and  reflect  the  new  volume  tructions- 


i2-l5) 


2-6 


NAVSWC  TR  91-718 
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FIGURE  2-1.  SLIC  INTERFACE  RECONSTRUCTION 
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CHAPTER  3 
NUMERICAL  RESULTS 

Calculations  of  shallow  depth  underwater  explosions  using  both  the  incompressible  and 
compressible  codes  are  presented  in  this  chapter.  Using  the  incompressible  code  a  parametric  study 
was  performed  on  the  charge  depth  which  ranged  from  5  to  20  feet.  The  validity  of  the  numerical 
results  was  checked  at  several  depths  where  the  calculations  were  repeated  on  a  computational  grid 
with  finer  resolution. 


INCOMPRESSIBLE  RUNS 


In  the  incompressible  code,  the  underwater  explosion  bubble  is  modeled  as  a  void  (zero  den¬ 
sity)  with  a  uniform  pressure  governed  by  the  adiabatic  gas  law 

PV'<  =  constant , 


w'here 


and 


P  is  the  bubble  pressure , 

V''  is  the  bubble  volume , 

y  =  1.3  is  the  ratio  of  specific  heats  of  the  bubble  gases . 


Initially,  the  bubble  is  assumed  to  be  at  rest  (that  is,  the  velocity  of  the  surrounding  water  is  zero) 
and  has  an  initial  radius  Rq  -  pressure  Pq.  Values  for  and  Pq  are  determined  using 


Po  =  PM-y) 


1-a^ 
1-fl  3(1-7) 


(3-1) 


(3-2) 


where 


J  -  13.1  is  the  charge  radius  constant  for  TNT, ^3 
N  -  0.023  is  the  charge  radius  ratio  for  TNT, 

W  =  100  is  the  charge  weight  in  pounds, 

P,,,  is  the  hydrostatic  pressure  at  the  charge  depth. 
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and 


a  = 


is  the  ratio  of  maxunmn  to  nunitrnun  bubble  radii. 


The  radius  ratio  N  is  an  empirically  determined  constant  such  that 


a 


1 


provided  that  is  measured  in  units  of  feet  of  water  (fw).  In  the  above,  is  the  “free  field" 
value  for  the  maximum  bubble  radius  in  the  case  of  a  bubble  in  an  infinite  medium  in  the  absence 
of  gravity.  (In  this  case  P ^  denotes  the  ambient  pressure  or  the  pressure  at  infinity.)  Given  a  value 
for  a,  and  an  ambient  pressure  P equation  (3-2)  is  derived  from  an  exact  integration  of  the  equa¬ 
tions  of  motion  for  a  spherically  symmetric  bubble  in  an  infinite  medium.  This  equation  is  often 
referred  to  as  the  Rayleigh  Piesset  equation.* 

Incompressible  calculations  were  performed  simulating  a  one  hundred  pound  charge  of  TNT  at 
depths  ranging  from  5  to  20  feet.  For  this  charge  it  follows  from  (3-1)  that  =  1  39851  feet.  The 
ambient  hydrostatic  pressure  is  determined  from 


where 

and 


Po.  =  c/  +  P, 

d  is  the  initial  charge  depth , 

=  33.9  feet  of  water  is  the  air  pressure . 


(3-4) 


Table  3-1  lists  the  values  of  the  initial  bubble  pressure,  free  field  ma.ximum  bubble  radius,  and 
scaled  charge  depth 


for  the  cases  studied  in  this  report. 
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TABLE  3-1.  SHALLOW  DEPTH  EXPLOSION  BUBBLE  INITIAL  DATA 


d  1 
(ft)  1 

^0  ! 
(fw)  j 

R  1 

(ft) 

5.0 

27401.4 

17.9 

0.279 

8.0 

27469.6 

17.5 

0.457 

10.0 

27513.5 

17.2  ‘ 

1 

0.580 

12.0 

27555.9 

17.0 

0.707 

14.0 

27597.2 

26.7 

0.838 

16.0 

27637.5 

16.5 

0.969 

18.0 

27676.7 

16.3  1 

1.104 

20.0 

27714.9  i 

16.1  ^ 

1.243 

Each  grid  used  for  the  calculations  consisted  of  a  region  of  uniform  cells  of  si.'e  Ar  =  A:  =  h 
in  the  region  r  <  27.0,  and  -36.0  <  z  <  72.0  of  a.xisymnieiric  (r.z)  space.  For  our  coarse  grid  runs 
this  region  was  divided  into  39x156  cells.  Cell  stretching  was  used  to  extend  the  boundaries  of  the 
computational  domain  to  z  =  -100,  z  =  140,  and  r  =  140.  The  grid  spacing  was  halved  in  each 
direction  for  the  fine  grid. 

Figure  3-1  displays  density  contours  and  velocity  vectors  at  six  different  times  for  the  coarse 
grid  run  at  charge  depth  d  =  5  feet.  In  this  figure,  ten  contours  are  shown,  with  five  concentrated 
betw'een  zero  and  0.1  Po  and  the  remaining  five  between  0.1  Po  and  Pq. 

In  the  first  frame  of  Figure  3-1  the  bubble  has  already  vented.  The  contours  above  the  cavity 
represent  regions  of  “spray”  w'here  the  density  is  below  the  value  (1-Ep)py=0.96p().  and  is  treated 
as  a  nonliquid  region  in  the  code.  Once  the  bubble  “vents”,  (a  bubble  cell  comes  in  contact  with 
an  air  cell)  the  pressure  of  the  cavity  is  given  the  value  of  the  ambient  air  pressure  instantaneously 
(see  Reference  9  for  details).  Therefore,  the  interaction  of  the  spray  and  the  high  pressure  explo¬ 
sive  gases  with  the  air  is  not  adequately  described  with  the  incompressible  model. 

The  second  frame  of  Figure  3-1  shows  the  cavity  at  its  maximum  depth  of  about  20  feet,  just 
before  it  begins  its  rebound.  The  initial  stages  of  the  rebound  is  seen  in  the  third  frame  (at  /  =  1 
second)  which  shows  the  cavity  refilling  from  the  bottom  up,  while  the  outer  fringes  of  the  cavity 
are  expanding  out  and  falling  away.  At  this  point  the  velocity  of  the  water  along  the  centerline 
attains  its  maximum  value  of  approximately  45  feet  per  second.  In  the  fourth  frame  a  tapered 
column  of  water  has  clearly  formed,  rising  approximately  10  feet  above  the  sea  level.  This  column 
continues  to  move  upw'ard  with  a  velocity  of  about  32  feet  per  second.  The  water  column  reaches 
its  maximum  height  of  25  feet  in  the  fifth  frame  (r  =;  2.4  seconds).  Here,  the  column  has  a  diameter 
of  approximately  34  feet  at  sea  level,  tapering  to  a  22-foot  diameter  at  a  10-fooi  elevation,  and  a 
15-foot  diameter  at  a  20-foot  elevation.  Finally,  the  sixth  frame  show's  the  column  beginning  to  fall 
downward,  due  to  the  influence  of  gravity. 
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Figure  3-2  displays  the  computed  evolution  of  tite  explosion  of  a  charge  at  a  dcpili  of  12  feci. 
The  first  frame  shows  an  early  stage  of  the  bubble  expansion.  I'he  surface  abo\c  the  bubble  has  a 
slight  bulge,  and  the  bubble  has  not  yet  vented  into  the  atmosphere.  The  contours  shown  inside  the 
bubble  are  a  numerical  tiriifact  due  to  remnants  of  cells  which  uere  iniiiaii)  only  paruall)  lilled 
with  water  (p  <  0.96po),  and  were  treated  as  nonliquid  cells  having  no  pressure  gradients,  and  there 
fore  no  velocity.  The  bubble  vents  at  r  =  0.2  second,  but  the  cavity  continues  to  increase  in  size, 
reaching  its  ma.ximum  depth  of  approximately  27  feet  in  the  second  frame  (7  =  0.4  second),  .o  thin 
layer  of  relatively  high  density  liquid  (not  adequately  resolved  on  tliis  grid)  persists  above  the  cav¬ 
ity,  even  until  the  third  frame  where  the  cavity  rebound  has  commenced  Similar  to  the  d  -  5  case, 
the  centerline  velocity  of  the  water  attains  its  maximum  value,  of  approximately  5.S  feet  per  second 
at  this  time.  The  water  column  has  risen  15  feet  above  sea  level,  and  is  moving  wuh  a  velocity  ol 
about  40  feet  per  second  in  the  fourth  frame  (r  =  1.5  seconds).  Tlie  tifih  frame  (i  ~  2.S  seconds') 
shows  the  water  column  at  its  maximum  height  of  42  feet.  This  column  has  a  base  diameter  of  45 
feet,  tapering  to  20  feet  at  a  20-foot  elevation,  and  approximately  maintaining  this  diameter  to  the 
top.  The  bulge  seen  at  a  30-foot  elevation  has  a  slightly  lower  density  than  water.  The  final  frame 
is  very  similar  to  the  previous  frame,  but  the  velocity  vectors  indicate  that  the  column  is  beginning 
to  fall  under  the  influence  of  gravity.  The  remnants  of  the  low  density  layer  can  still  be  seen  above 
the  column. 

The  d  =  12  case  was  repeated  on  the  fine  grid  as  a  numerical  validity  check  The  fine  grid 
results  shown  in  Figure  3-3  reveal  more  detail  in  the  “domed”  region  above  the  bubble  at  early 
times,  a  greater  resolution  of  the  outer  thin  layer  surrounding  the  water  column  at  later  times,  and  a 
thin  jet  of  water  and  spray  ejected  upward  along  the  centerline  axis.  The  water  column  shape, 
width,  and  height  predictions  are  very  similar  for  the  two  computations. 

A  more  quantitative  comparison  between  the  coarse  and  fine  grid  compulations  of  the  d  =  12 
case  is  presented  in  Figure  3-4.  In  this  figure  the  column  height  is  defined  to  be  the  location  of  the 
lowest  contour  of  p  =  0.96po  crossing  the  axis,  the  uatcr  height  is  the  highest  contour  of 
p  =  0.96po,  and  the  spray  height  is  the  highest  contour  of  p  =  0.1  Py.  The  column  heights  (which  at 
early  times  represent  the  location  of  the  bottom  of  the  bubble)  are  in  close  agreement  throughout  the 
entire  time  interval.  The  discrepancies  in  the  water  heights  indicate  that  the  coarse  grid  has 
insufficient  resolution  for  the  top  of  the  domed  region  above  the  bubble  and  subsequently  above  the 
w'ater  column.  In  the  coarse  grid  computation  no  water  cells  remain  above  the  water  column  after 
approximately  t  =  1.4  seconds.  The  jagged  graphs  of  both  the  water  height  and  spray  height  from 
the  coarse  grid  calculation  are  also  due  to  insufficient  resolution  of  the  dome  region.  The  graph  of 
the  maximum  spray  height  for  the  fine  grid  follows  a  nearly  parabolic  trajectory  as  it  is  decelerated 
by  gravity,  reaching  a  maximum  height  of  over  80  feet. 

The  evolution  of  an  explosion  with  an  initial  charge  depth  of  16  feet  is  displayed  in  Figure  3-5. 
In  the  second  frame  the  bubble  attains  its  maximum  size.  A  significant  “bump”  has  fomied  above 
the  bubble,  having  an  average  density  of  approximately  0  Spy.  The  code  treats  this  bump  as  a  uni¬ 
form  pressure  nonliquid  region  being  decelerated  by  gravity.  This  bump  is  caused  by  the  initial 
upward  acceleration  due  to  the  high  pressure  bubble,  followed  by  a  deceleration  once  the  bubble 
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pressure  falls  below  the  atmospheric  pressure.  In  the  third  frame  the  water  jet  has  passed  thiough 
the  bubble  and  is  impacting  a  thin  layer  of  water  which  was  the  top  surfttce  of  the  bubble  At  this 
point  the  bubble  is  an  annular  region.  The  founh  frame  shows  the  svater  jet  rising  as  a  column  of 
water  about  twenty  feet  above  the  surface.  The  water  column  is  being  impeded  by  the  relatively 
high  density  spray  in  the  bump  which  is  now  falling  downward.  The  interaction  of  the  water 
column  with  the  spray  causes  a  radially  outward  movement  of  spray  as  depicted  in  tire  liftii  frame. 
Here,  the  water  column  has  reached  its  ma.ximurn  height  of  about  40  feet,  but  its  diameter  has 
diminished  to  only  10  feet.  However,  a  significant  amount  of  spray,  having  a  density  of  about  0.2py 
has  spread  out  in  a  region  having  an  80-foot  diameter.  Finally,  the  sixth  frame  shows  the  column 
and  spray  falling  downward  due  to  gravity. 

The  d  =  \6  case  was  also  repeated  on  the  fine  grid.  These  results  are  shown  in  Figure  3-6, 
which  are  qualitatively  similar  to  the  coarse  grid  results  (Figure  3-5)  for  the  first  four  frames,  but 
reveal  a  much  higher  and  thinner  water  column  in  the  final  two  frames.  F’unhennore,  tins  coiiimn  is 
surrounded  by  a  complex  pattern  of  relatively  den.se  layers  having  a  maximum  diameter  of  about  50 
feet,  and  rising  to  an  elevation  of  about  60  feet.  The  complexity  of  the  flow-  field  at  the  later  times 
is  due  to  the  interactions  of  the  primary  water  jet  with  the  top  surface  of  the  bubble. 

Quantitative  differences  between  the  coarse  and  fine  grid  results  for  the  d  =  \6  case  are  shown 
in  Figure  3-7.  The  almost  discontinuous  increases  in  the  column  heights  (panicularly  evident  for 
the  fine  grid  run)  are  caused  when  the  water  jet  impacts  the  top  pan  of  the  bubble,  which  has 
formed  a  downward  jet.  This  downward  jet  can  be  seen  in  the  second  frame  of  Figure  3-6.  This 
phenomenon  has  also  been  reponed  experimentally  by  Blake  and  Gibson and  computationally 
using  a  boundary  integral  method  by  Blake,  Taib  and  Doherty.'^  Just  after  the  impact  of  the  top  and 
bottom  surfaces  of  the  bubble,  the  column  height  is  abruptly  changed  from  the  height  of  the  lower 
jet  to  the  top  of  the  water  level  along  the  axis.  The  relatively  large  discrepancies  in  the  spray  and 
column  heights  at  later  times  suggest  that  further  grid  refinement  studies  should  be  performed  due  to 
the  complex  nature  of  this  problem. 

Figures  3-8  and  3-9  show  the  time  histories  of  the  water  column  heights  for  the  coarse  grid 
runs  at  depths  betw-een  5  and  14  feet,  and  15  and  20  feet,  respectively.  The  steepness  of  the  cur\es 
in  the  deeper  cases  is  once  again  due  to  the  lower  jet  impacting  the  water  at  the  top  of  the  bubble. 
In  the  shallow'  cases  the  water  at  the  top  of  the  bubble  has  been  pushed  away  due  to  the  initial 
acceleration  of  the  explosion.  These  figures  indicate  that  a  depth  of  14  feet  is  optimal  for  attaining 
the  highest  water  column.  However,  while  the  shallow  depth  cases  are  probably  being  predicted 
accurately  (as  evidenced  by  the  small  discrepancies  in  the  coarse  and  fine  grid  runs  at  12  feet),  the 
errors  in  the  deeper  cases  are  greater  and  higher  column  heights  than  that  predicted  from  the  coarse 
grid  runs  can  be  expected  (cf.  Figure  3-7),  However,  it  is  probably  safe  to  conclude  that  the 
optimal  depth  for  attaining  the  maximum  column  height  is  between  12  and  16  feet  for  a  lOO-pound 
charge  of  TNT. 

The  time  histories  of  the  computed  coarse  grid  spray  heights  are  displayed  in  Figures  3-10  and 
3-]  I.  The  abrupt  drop  in  the  maximum  spray  height  values  at  the  shallow  charge  depths  (less  than 
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twelve  feet)  are  primarily  numerical  artifacts  due  to  the  grid  streichiiig  abu\e  '^2  feet  A*-  ihc  '>pra) 
moves  upward  into  the  larger  cells,  the  average  cell  densities  decrease  until  die;.'  tall  bciuu  O.lp,, 
Table  3-2  lists  extrapolated  spray  heights  obtained  from  tlie  spray  heights  before  the  abrupt  drop  oft 
The  extrapolation  is  based  on  the  assumption  that  tlie  spray  is  being  decelerated  b>  gravity  Since 
wind  resistance  is  neglected  in  these  extrapolations,  vve  expect  these  values  are  l.irger  I'n.m  those 
which  would  result  from  actual  experiments. 


TABLE  3-2.  EXTRAPOLATED  SHALLOW  CHARGE  DEPTH  SPRAY  HEIGHTS 


. . 

.-Xpproximale 

Initial 

X'elocity 

(ft 

Maximum 

Spray 

Heistht 

sec ) 

Time  i 

: 

Maximum  i 

(ft)  i 

5 

300 

1-100 

9.4 

8 

180 

500 

5.6 

10 

120 

220 

3.7 

12 

80 

100 

Experimental  data  tor  3()0-pound  charges  of  I’.X'T  itas  been  collected  by  Young.*  Since  the  bubble 
scales  with  the  cube  root  of  the  charge  weight  (cf.  3-1)  our  results  can  be  compared  to  those  in 
Reference  2  by  scaling  the  depth  by  3'  \  Except  for  the  extrapolated  heights  for  tlie  shallow  causes 
u=5  and  u-h.  ihe  results  in  Ttible  3-2  and  Figures  3  10  and  3-11  greatly  under  estimate  the  max¬ 
imum  plume  heights  reported  by  Young. However,  we  are  unaware  of  publications  reporting  the 
density  inside  of  experimentally  produced  plumes  Jt  i^  possible  that  they  are  compri.scd  of  a  very 
low  density  mix  of  steam  and  explosion  gases,  which  are  not  modeled  using  the  incompressible 
code. 

COMPRESSIBLE  RUNS 

Two  computations  were  carried  out  using  the  compressible  code.  The  first  simulated  a  3()() 
pound  TNT  explosion  at  a  depth  of  5  feet.  Here  a  perfect  gas  equtiiion  of  state  was  used  to  model 
the  air/e.xplosive  products  with  y  =  1.4.  This  calculation  was  carried  out  in  two  phases 

1.  A  1-D  spherical  calculation  tracked  the  initial  expansion  of  the  explosion  bubble  to  the  point 
where  the  explosive  shock  neared  the  water  surface. 

2.  A  2-D  axisymmetric  calculation  continued  the  calculation  to  later  times 
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The  i-D  calculation  was  initiated  from  a  constant  pressure  bubble  at  rest  uitii  the  t'ullowine  proper¬ 
ties; 


Pq  =  7.6742(10'^')  dynes icni',  p  =  1.63  yrnicm^,  R^-,  ~  26  84  crn 

The  results  of  this  computation  are  shown  in  Figures  3-12  and  3-13.  Figure  3-12  illusirates  the 
velocities  in  both  air  and  water  throughout  the  simulation.  The  heavy  line  running  through  these 
figures  is  the  air-water  interface.  The  complex  interface  geometry  is  more  clearly  visible  in  Figure 
3-13  which  provides  a  density  contour  plot.  The  description  of  the  explosion  provided  by  the  simu¬ 
lation  can  be  divided  into  four  phases; 

1.  Bubble  Venting.  The  bubble  vents  within  the  first  10  to  20  milliseconds. 

2.  Cavity  formation.  A  deep  air  cavity  is  carved  out  in  the  water  by  the  explosive.  This  cavity 
reaches  its  maximum  depth  at  about  0.5  second. 

3.  Cavity  collapse.  Cavity  venting  leads  to  a  gas  pressure  which  is  too  low  to  sustain  the  cavity. 
Water  rushes  into  the  cavity  to  re-establish  to  original  water  pressure. 

4.  Rebound.  Water  rushing  into  the  cavity  impacts  at  the  cavity  center  producing  a  plume,  which 
in  this  case  rises  to  a  height  of  4  meters 

Absent  from  the  above  description  is  the  interaction  of  the  original  explosion  shock  with  the  sur¬ 
face.  In  this  case,  the  time  difference  between  venting  and  the  arrival  of  the  shock  is  very'  small, 
and  the  shock  does  not  appear  to  have  a  significant  effect  on  the  plume. 

The  second  case  considered  consisted  of  100  pounds  of  TNT  exploded  at  a  depth  of  6.5  feet  in 
water  which  was  13.1  feet  deep.  For  this  case,  the  explosive  producis/air  were  modeled  by  the 
Jones-Wilkins-Lee  (JWL)  equation  of  state,  which  has  the  fonn: 

p  =/i(v)  +/2(v)  +  — 


(3-5) 


For  TNT  the  constants  are  defined  as: 

■U]  =  3.760(10’^)  dynes /cni^ 

02  =  3.273(10'^^)  dynes /cm" 

C,  =  -.0443492  an^/gm 

C^  -  -.19373  cm^/gm  (3-6) 

X\  =  -6.7645  gm/crn^ 

X2  -  -1.5485  gm/cm^ 

CO  =  .30. 
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The  calculation  was  initiated  using  a  one  dimensional  detonation  model'^’  winch  procidcd  a  saii: 
unifomi  pressure,  velocity  and  density  profile  at  the  instant  at  which  tlie  dcionatatn  iruni,  uinefi  a. is 
initiated  at  the  center  of  tlie  charge,  reached  the  charge  suiface.  Tlie  radius  of  the  initial  vliargc  w.t- 
18.8  cm  .As  in  the  previous  calculation,  this  case  was  continued  using  a  1  1)  soiunon  iintd  me 
shock  was  close  to  the  surface.  The  calculation  was  then  restarted  using  2  D.  asMiupctric  c>.Hir<li 
nates. 

The  results  from  the  second  calculation  are  shown  in  Figures  .3-14  and  .v  l5  Figure  .F!4  dUis 
trates  the  velocity  vectors  for  both  air  and  water,  with  the  interface  siiown  as  a  solid  line  Dcnsits 
contours  arc  given  in  Figure  3-15  and  clearly  indicate  the  evolution  of  tlie  water-gas  interface,  Tiic 
explosion  phases  shown  in  Figures  3-14  and  3-15  clearly  parallel  those  of  itie  tirsi  case  (Figures  3- 
12  and  3-13).  The  venting,  cavity  formation,  and  cavity  collapse  stages  arc  a'!  esidcnt.  h.owcser,  a 
strong  rebound  is  not.  This  may  be  a  consequence  of  the  large  region  of  air- water  mixture  which 
appears  along  the  center  line.  Presumably,  the  entrained  air  imparts  a  compressibility  to  the  mnxiure 
which  absorbs  the  energy  of  cavity  collapse,  thus  preventing  the  formation  of  a  w.iicr  column. 
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F1(;LRE  3-4.  TIME  HISTORIES  OF  THE  MAXIMUM  WATER  AND  SPRAY 
HEIOHTS  FOR  THE  TWELVE-FOOT  DEPTH  COMPUTATIONS 
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FIGURE  3-7.  TIME  HISTORIES  OF  THE  MAXIMUM  COl^UMN  AND  SPRAY 
HEIGHTS  FOR  THE  SIXTEEN-FOOT  DEPTH  COMPUTATIONS 
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FIGURE  3-8.  TIME  HISTORIES  OF  THE  MAXIMUM  COLUMN  HEIGHTS 

BASED  ON  THE  COARSE  GRID  COMPUTATIONS  AT  CHARGE 
DEPTHS  BETWEEN  FIVE  AND  FOURTEEN  FEET 
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FIGURE  3-9.  TIME  HISTORIES  OF  THE  MAXIMUM  COLUMN  HEIGHTS 

BASED  ON  THE  COARSE  GRID  COMPUTATIONS  AT  CHARGE 
DEPTHS  BETWEEN  FIFTEEN  AND  TWENTY  FEET 
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FIGURE  3-10.  TIME  HISTORIES  OF  THE  MAXIMUM  SPRAY  HEIGHTS 

BASED  ON  THE  COARSE  GRID  COMPUTATIONS  AT  CHARGE 
DEPTHS  BETWEEN  FIVE  AND  FOURTEEN  FEET 
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FIGURE  3-11.  TIME  HISTORIES  OF  THE  MAXIMUM  SPRAY  HEIGHTS 

BASED  ON  THE  COARSE  GRID  COMPUTATIONS  AT  CHARGE 
DEPTHS  BETWEEN  FIFTEEN  AND  TWENTY  FEET 
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FIGURE  3-12.  (CONT.)  VELOCITY  VECTORS  FROM  THE  COMPRESSIBLE  COMPUTATIO.N 
OF  A  300-POUND  CHARGE  OF  TNT  AT  A  DEPTH  OF  FIVE  FEET 
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FIGURE  3-13.  DENSITY  CONTOURS  FROM  THE  COMPRESSIBLE  COMPUTATION 
OF  A  300-POUND  CHARGE  OF  TNT  AT  A  DEPTH  OF  FIVE  FEET 
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FIGURE  3-14.  VELOCITY  VECTORS  FROM  THE  COMPRE: 
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[CONT.)  VEl.OCi  rY  VECTOR.S  FROM  THE  COMPRES.SIBi.E  COMI’IJT.V  I'lON  OF 
lOO-POUND  CIIAROE  OF  TN!'  IN  SVIAl.EOVV  WATER  AI  A  DEPm  OF  6.5  FEET 
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FIGURE  3-15.  DENSITY  CONTOURS  FROM  THE  COMPRESSIBLE  COMPUTATION 
OF  A  lOO-POUND  CHARGE  OF  TNT  IN  SHALLOW  WATER 
AT  A  DEPTH  OF  6.5  FEET 
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CHAPTER  4 
CONCLUSIONS 

In  this  paper  two  different  methods  were  described  and  used  to  predict  the  dynamics  of  shal¬ 
low  depth  explosion  plumes.  The  incompressible  code  has  the  advantage  of  computational 
efficiency  primarily  because  it  selects  time  steps  based  on  the  velocity  of  the  water  rather  than  its 
much  larger  sound  speed.  This  code  also  has  the  advantage  of  being  validated  on  predicting  under¬ 
water  bubble  dynamics  at  greater  depths.^’*’^  The  compressible  code  has  the  advantage  of  a  better 
physical  model  in  which  the  effects  of  the  initial  shock  and  the  subsequent  interaction  of  air  and 
water  can  be  predicted. 

At  early  times  after  a  shallow  depth  underwater  explosion,  the  water  above  the  bubble  is 
pushed  radially  outward,  forming  a  thin  water  dome  rising  above  the  initial  water  level.  A  critical 
event  which  may  occur  during  this  time  is  the  venting  of  the  bubble  through  the  thin  water  dome 
into  the  atmosphere.  Both  the  incompressible  and  compressible  code  calculations  presented  in  this 
report  predict  that  the  bubble  vents  into  the  atmosphere  before  the  occurrence  of  the  first  maximum 
bubble  volume.  The  venting  of  the  bubble  will  be  predicted  numerically  by  the  incompressible  code 
whenever  the  density  of  a  cell  containing  the  thin  dome  falls  below  (]-ep)po  (see  (2-2)).  Thus,  the 
numericai  prediction  of  venting  is  a  function  of  both  the  grid  size  and  the  value  ep  used  to  deter¬ 
mine  the  cutoff  between  the  liquid  and  nonliquid  cells.  Preliminary  numerical  studies  have  shown 
that  the  selection  of  Ep  can  effect  the  prediction  of  venting,  and  have  a  great  effect  on  the  subse¬ 
quent  dynamics  of  the  bubble,  cavity  and  plume.  These  preliminary  computations  are  in  better 
agreement  with  the  experimental  results  of  Blake  and  Gibson, which  show  that  spark  generated 
bubbles  will  not  vent  at  values  of  00.56.  Shallow  depth  explosion  bubbles  which  do  not  vent  form 
a  thin  water  jet  above  the  dome  instead  of  the  large  cavity  which  opens  up  for  a  venting  bubble. 
The  incompressible  simulations  presented  here  predict  venting  at  early  times  for  values  of  c  as  large 
as  1.  Because  of  this,  the  reliability  of  the  calculations  could  be  greatly  enhanced  by  empirically 
selecting  appropriate  values  for  Ep  for  the  incompressible  code  to  better  match  the  predictions  of 
venting.  Subsequent  research  in  developing  improved  predictive  capabilities  will  be  focused  on  this 
aspect  of  the  calculations. 
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